addpath './KernelGradient'
addpath './BoundaryIntegral'
addpath './util'
addpath './ElTopo/'

cx = [-1 0 1 1 -1];
cy = [1 0.8 1 -1 -1];

% bx = [-1 1 1 -1];
% by = [1 2 0 -1];

[x, y] = generateBoundaryNode(cx, cy, 30);
% [dx, dy] = generateBoundaryNode(bx, by, 100);

domain2D = createDomain(x, y);

dim = size(domain2D.x, 2);

% b = zeros(2*dim, 1);
% b(1:2:2*dim) = dx' - x';
% b(2:2:2*dim) = dy' - y';

b = assemblyBoundaryVector(domain2D);

res = 50;
[rx, ry] = meshgrid(0:res,0:res);
rx = rx / res * 1.86 - 0.93;
ry = ry / res * 1.86 - 0.93;
rx = reshape(rx, (res+1)*(res+1), 1);
ry = reshape(ry, (res+1)*(res+1), 1);

% rx = [-0.5; -0.5; 0.5; 0.5];
% ry = [-0.5; 0.5; 0.5; -0.5];

[A, Z] = assemblyMatrix(domain2D);

% b = assemblyBoundaryVector(domain2D);

theta = A \ b;

u = evaluateInteriorPoint(rx, ry, domain2D, theta);

F = evaluateInteriorDeformationGradient(rx, ry, domain2D, theta);

plotDeformation(domain2D, rx, ry, b(1:2:2*dim), b(2:2:2*dim), u(:,1), u(:,2));

% dx = rand(1, 1);
% dy = rand(1, 1);
% 
% rc = rand(1000*9, 1) * 100;
% strain = reshape(rc, 3, 3, 1000);
% 
% plotDeformation(domain2D, rx, ry, zeros(dim, 1), zeros(dim, 1), dx, dy);